Does modernization lead to secularization in societies?
Let’s unpack some definitions:
Modernization refers, broadly, to economic development and institutions with some relation to economic development. Per capita GDP (PPP) is an obvious direct measurement of this, but it is also incomplete: a country may have a high average due to the presence of extractable natural resources that often lead to extremely skewed wealth distributions and, as such, averages that aren’t necessarily reflective of overall development. Life expectancy, safety standards, etc. also reflect modernization.
Secularization refers to a general shift away from religious practice and religiously-based values. This can be a bit nebulous as different faiths may emphasize different values, and is particularly problematic when comparing countries with different historical faith supersets (e.g. Islamic vs. Christian vs. Buddhist). For this analysis, I’ll include only countries in Western and Eastern Europe, North and Latin America, and Australia/New Zealand. While it would be beneficial to include sub-Saharan African countries, insufficient data exists in a number of public opinion-related questions, and as such, inference on them would rely heavily on imputation.
The basic idea of the modernization hypothesis is that as countries become wealthier, their citizens become more secular/less religious. The nature of this mechanism (by mechanism, we could essentially think of why it is that a maps to b) isn’t exactly universally agreed upon. For example, an economic theory might say that religion is a form of social insurance (churches often provide social services, and faith may in and of itself also act as an intangible insurance–a future payout for enduring hardships of today) and that when the state and industries within the state provide this additional insurance through better public and private goods, demand for the good diminishes. For our purposes, I’ll concentrate less on whether that explanation or another makes more sense and concentrate on the empirics–all hail data science!
For this endeavour I collected data from a number of sources–some on wikipedia, others in reports.
Constructing measurement models is a hobby of mine.
First, let’s scale the data. Scaling is important because variables measured on very different scales will lead to spurious results, due to the larger weight exerted by variables captured in numerically larger values. There are several scaling methods to choose from, but I’ll go with a standard conversion to a standard normal variable (e.g. distributed mean 0, standard deviation 1), through subtracting the mean and dividing by the standard deviation. Note here that I coded variables in such a manner as to presumably achieve positive correlations as opposed to a mix of positive and negative correlations–so sometimes ranges are [0, -inf) as opposed to [0, +inf). This just makes a few steps easier down the road.
Note that you’ll see a subsetting of the data: I initially collected data for more regions than I ultimately opted to include. You’ll also see an imputation step due to limited availability for some measurements. The countries with the largest missingness problems (e.g. Pacific Islands) are excluded from eventual analysis. This leaves 64 countries with 53 total
##
## iter imp variable
## 1 1 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 1 2 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 1 3 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 1 4 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 1 5 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 2 1 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 2 2 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 2 3 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 2 4 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 2 5 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 3 1 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 3 2 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 3 3 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 3 4 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 3 5 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 4 1 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 4 2 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 4 3 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 4 4 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 4 5 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 5 1 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 5 2 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 5 3 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 5 4 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
## 5 5 Religion.Important PubOp Church Suic Water COL Internet Safety Specialists DocSal DocRelSal Irreligion Informal.Payments Power.Outages Prosperity.Index Bribes Mercer Innovation RuleOfLaw Satisfaction CapRent CountryAptCost SocialCap AdulteryUnac HomUn AbortionOp PremSexUnac GAI
To get an early hint of the underlying dimensionality of the data, let’s perform a singular value decomposition and visualize the singular values. In general this will resemble an elbow. If there was one latent variable reasonably capturing all of these different measurements, then there would be one large value followed by k-1 small ones, where k is the total number of variables in the data.
svd(df[,c(3:k)])$d %>% plot(pch=16, main="Singular Value Decomposition of Modernization Variables", xlab="Singular Value", ylab="Magnitude of Singular Value")
As another way of visualizing, let’s do a correlation plot.
corrplot(cor(df[,c(3:k)]))
Ok, now for some fun stuff. And by fun I mean “math that R will do for me.”
When I first envisioned this endeavor, I preferred creating separate scales for governmental development and economic development. However, I quickly realized that separating the two did not make sense. The development of governmental institutions (reflected in metrics such as Rule of Law and Judicial Independence from the Index of Economic Freedom, and World Bank data such as informal payments being made) correlate so strongly with measurements of economic/social development (life expectancy, GDP Per Capita) that concentrating on the two separately would introduce redundancy. It would be redundantly redundant (wink wink).
As such, I examined a reliability measurement (Cronbach’s Alpha) of the selected modernization variables. And as you’ll see below, the alpha is pretty good, well above the informal threshold for what constitutes a reasonable scale. This isn’t exactly the most surprising thing though–is it a shocker that countries with more bribes have less governmental integrity? Some of these are essentially different organizations’ measurements of the same thing. Luckily, while correlated, they diverge sufficiently so as to not introduce mathematical redundancy, which would prevent matrices from being able to be inverted along with other such horrors.
##
## Reliability analysis
## Call: psych::alpha(x = df %>% select(Life.Expectancy, GDP.Per.Capita,
## Internet, Prosperity.Index, Innovation, RuleOfLaw, Integrity,
## Judicial, Pollution, Power.Outages, AirPollDeaths, Informal.Payments,
## Bribes, Fragility, SocialCap, HospitalBeds))
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.96 0.96 0.98 0.62 26 0.0065 3.2e-17 0.8 0.66
##
## lower alpha upper 95% confidence boundaries
## 0.95 0.96 0.98
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r
## Life.Expectancy 0.96 0.96 0.98 0.62 24 0.0070 0.055
## GDP.Per.Capita 0.96 0.96 0.98 0.61 23 0.0073 0.055
## Internet 0.96 0.96 0.98 0.63 26 0.0066 0.057
## Prosperity.Index 0.96 0.96 0.98 0.60 22 0.0075 0.051
## Innovation 0.96 0.96 0.98 0.61 23 0.0073 0.054
## RuleOfLaw 0.96 0.96 0.98 0.60 23 0.0075 0.052
## Integrity 0.96 0.96 0.98 0.60 23 0.0075 0.052
## Judicial 0.96 0.96 0.98 0.61 23 0.0072 0.054
## Pollution 0.96 0.96 0.98 0.63 25 0.0068 0.055
## Power.Outages 0.96 0.96 0.98 0.64 27 0.0063 0.056
## AirPollDeaths 0.96 0.96 0.98 0.63 25 0.0067 0.053
## Informal.Payments 0.96 0.96 0.98 0.64 27 0.0064 0.052
## Bribes 0.96 0.96 0.98 0.61 24 0.0072 0.056
## Fragility 0.96 0.96 0.98 0.60 23 0.0074 0.053
## SocialCap 0.96 0.96 0.98 0.61 24 0.0071 0.057
## HospitalBeds 0.97 0.97 0.98 0.68 32 0.0055 0.025
## med.r
## Life.Expectancy 0.65
## GDP.Per.Capita 0.65
## Internet 0.69
## Prosperity.Index 0.65
## Innovation 0.65
## RuleOfLaw 0.65
## Integrity 0.65
## Judicial 0.65
## Pollution 0.66
## Power.Outages 0.69
## AirPollDeaths 0.67
## Informal.Payments 0.69
## Bribes 0.66
## Fragility 0.65
## SocialCap 0.66
## HospitalBeds 0.69
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## Life.Expectancy 64 0.84 0.84 0.84 0.81 -5.8e-17 1
## GDP.Per.Capita 64 0.90 0.90 0.90 0.88 5.7e-19 1
## Internet 64 0.71 0.71 0.69 0.67 -1.4e-16 1
## Prosperity.Index 64 0.97 0.97 0.98 0.97 6.2e-16 1
## Innovation 64 0.91 0.91 0.91 0.89 2.8e-16 1
## RuleOfLaw 64 0.95 0.95 0.96 0.94 6.3e-17 1
## Integrity 64 0.95 0.95 0.96 0.94 1.3e-16 1
## Judicial 64 0.89 0.89 0.89 0.88 -1.4e-16 1
## Pollution 64 0.76 0.76 0.75 0.73 -1.4e-16 1
## Power.Outages 64 0.64 0.64 0.62 0.59 9.0e-18 1
## AirPollDeaths 64 0.73 0.73 0.73 0.69 -8.0e-17 1
## Informal.Payments 64 0.62 0.62 0.59 0.57 8.2e-17 1
## Bribes 64 0.87 0.87 0.86 0.85 -1.6e-17 1
## Fragility 64 0.94 0.94 0.94 0.93 1.9e-17 1
## SocialCap 64 0.86 0.86 0.85 0.84 -2.7e-16 1
## HospitalBeds 64 0.29 0.29 0.25 0.21 8.0e-17 1
Ok, convinced of the merits of my lovely scale? I sure hope so. Onwards! Let’s actually convert this to a single measurement. We’ll use what is probably the most common method around–exploratory factor analysis.
df5 <- df %>% select(
Life.Expectancy, GDP.Per.Capita,Internet,
Prosperity.Index,Innovation,RuleOfLaw, Integrity, Judicial,
Power.Outages,AirPollDeaths,Bribes, Fragility,
SocialCap, Mercer, Nurses, HospitalBeds
)
row.names(df5) <- df$Country
Development <- factanal(df5, factors=1, scores='regression')
Development
##
## Call:
## factanal(x = df5, factors = 1, scores = "regression")
##
## Uniquenesses:
## Life.Expectancy GDP.Per.Capita Internet Prosperity.Index
## 0.338 0.181 0.536 0.021
## Innovation RuleOfLaw Integrity Judicial
## 0.154 0.059 0.069 0.196
## Power.Outages AirPollDeaths Bribes Fragility
## 0.696 0.504 0.253 0.092
## SocialCap Mercer Nurses HospitalBeds
## 0.293 0.139 0.361 0.948
##
## Loadings:
## Factor1
## Life.Expectancy 0.814
## GDP.Per.Capita 0.905
## Internet 0.681
## Prosperity.Index 0.989
## Innovation 0.920
## RuleOfLaw 0.970
## Integrity 0.965
## Judicial 0.897
## Power.Outages 0.552
## AirPollDeaths 0.704
## Bribes 0.864
## Fragility 0.953
## SocialCap 0.841
## Mercer 0.928
## Nurses 0.800
## HospitalBeds 0.228
##
## Factor1
## SS loadings 11.158
## Proportion Var 0.697
##
## Test of the hypothesis that 1 factor is sufficient.
## The chi square statistic is 305.96 on 104 degrees of freedom.
## The p-value is 9.15e-22
Development$scores
## Factor1
## Albania -0.88064231
## Argentina -0.42545636
## Armenia -0.67876639
## Australia 1.22086443
## Austria 1.18512426
## Belarus -0.77869809
## Belgium 1.03196382
## Belize -1.04932109
## Bolivia -1.47655323
## Bosnia and Herzegovina -0.86088639
## Brazil -0.69985438
## Bulgaria -0.35213578
## Canada 1.24703984
## Chile 0.25878189
## Colombia -0.85041962
## Costa Rica 0.06906005
## Croatia -0.09581506
## Cyprus 0.21747533
## Czech Republic 0.53775705
## Denmark 1.62653447
## Ecuador -1.00036184
## El Salvador -1.12527114
## Estonia 0.91783286
## Finland 1.30588674
## France 0.94141051
## Georgia -0.50080105
## Germany 1.31351066
## Greece -0.05057861
## Guatemala -1.36071022
## Honduras -1.47863551
## Hungary -0.16110829
## Iceland 1.34551979
## Ireland 1.18637562
## Italy 0.42929085
## Latvia 0.14359776
## Lithuania 0.32928796
## Mexico -0.83548300
## Moldova -1.06405667
## Montenegro -0.37329942
## Netherlands 1.42089037
## New Zealand 1.38094200
## Nicaragua -1.51419947
## North Macedonia -0.72235418
## Norway 1.64040313
## Panama -0.45815580
## Paraguay -1.02238349
## Peru -0.81031461
## Poland 0.21142643
## Portugal 0.67466322
## Romania -0.11687965
## Russia -0.88867547
## Serbia -0.56701237
## Slovakia 0.13489740
## Slovenia 0.63539963
## Spain 0.74323570
## Suriname -0.99024716
## Sweden 1.52899112
## Switzerland 1.57263984
## Ukraine -0.99415617
## United Kingdom 1.13670182
## United States 0.95362966
## Uruguay 0.22013860
## Venezuela -2.22790680
## Guyana -1.15013315
I next create a scale for secularism. Well, I call it faith, but that’s just because it comes back decreasing in secularism and multiplying by -1 is arbitrarily out of the question. But basically, more secularization = less faith and vice versa. Below I select variables, calculate Cronbach’s Alpha to assess scale reliability, and view the results of the singular value decomposition of the matrix of values as an additional robustness check, because we all need a little more SVD in our lives.
psych::alpha(df %>% select(Religion.Important,PubOp,Church,Abor,Irreligion,
AdulteryUnac,HomUn,AbortionOp,PremSexUnac, GAI))
## Number of categories should be increased in order to count frequencies.
##
## Reliability analysis
## Call: psych::alpha(x = df %>% select(Religion.Important, PubOp, Church,
## Abor, Irreligion, AdulteryUnac, HomUn, AbortionOp, PremSexUnac,
## GAI))
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.91 0.91 0.94 0.49 9.8 0.018 -2.7e-17 0.74 0.51
##
## lower alpha upper 95% confidence boundaries
## 0.87 0.91 0.94
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r
## Religion.Important 0.89 0.89 0.93 0.48 8.4 0.020 0.035
## PubOp 0.90 0.90 0.93 0.49 8.7 0.019 0.033
## Church 0.90 0.90 0.93 0.49 8.6 0.020 0.038
## Abor 0.91 0.91 0.93 0.52 9.6 0.018 0.034
## Irreligion 0.89 0.89 0.92 0.48 8.1 0.021 0.033
## AdulteryUnac 0.91 0.91 0.95 0.54 10.7 0.016 0.027
## HomUn 0.89 0.89 0.93 0.47 8.1 0.021 0.036
## AbortionOp 0.89 0.89 0.93 0.47 7.8 0.022 0.034
## PremSexUnac 0.90 0.90 0.94 0.49 8.8 0.020 0.040
## GAI 0.90 0.90 0.93 0.51 9.2 0.018 0.025
## med.r
## Religion.Important 0.51
## PubOp 0.51
## Church 0.48
## Abor 0.55
## Irreligion 0.46
## AdulteryUnac 0.57
## HomUn 0.46
## AbortionOp 0.45
## PremSexUnac 0.51
## GAI 0.51
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## Religion.Important 64 0.79 0.79 0.79 0.73 -6.5e-19 1
## PubOp 64 0.75 0.75 0.74 0.68 1.2e-18 1
## Church 64 0.76 0.76 0.74 0.70 -8.8e-17 1
## Abor 64 0.63 0.63 0.60 0.54 1.9e-17 1
## Irreligion 64 0.83 0.83 0.83 0.78 5.0e-18 1
## AdulteryUnac 64 0.49 0.49 0.41 0.38 -1.3e-17 1
## HomUn 64 0.83 0.83 0.81 0.78 -9.8e-18 1
## AbortionOp 64 0.87 0.87 0.87 0.84 2.6e-17 1
## PremSexUnac 64 0.74 0.74 0.70 0.67 -6.9e-18 1
## GAI 64 0.68 0.68 0.67 0.59 -2.1e-16 1
svd(df %>% select(Religion.Important,PubOp,Church,Abor,Irreligion,
AdulteryUnac,HomUn,AbortionOp,PremSexUnac))$d %>%
plot(pch=16, main="Results of Singular Value Decomposition of Secularization Variables", xlab="Singular Value", ylab="Magnitude of Singular Value")
Ok, Cronbach’s Alpha looks ok (a bit less impressive than the first one but hey it has a lot of human questions in it and humans aren’t necessarily so consistent in their answers), and demonstrated with the SVD that one latent variable captures a good bit of the action.
With that being done, let’s rinse and repeat with some more factor analysis.
Faith <-df %>% select(Religion.Important,PubOp,Church,Abor,Irreligion,
AdulteryUnac,HomUn,AbortionOp,PremSexUnac)
row.names(Faith) <- df$Country
Faith <- factanal(Faith, factors=1, scores='regression')
Faith
##
## Call:
## factanal(x = Faith, factors = 1, scores = "regression")
##
## Uniquenesses:
## Religion.Important PubOp Church Abor
## 0.300 0.561 0.511 0.685
## Irreligion AdulteryUnac HomUn AbortionOp
## 0.244 0.857 0.426 0.210
## PremSexUnac
## 0.524
##
## Loadings:
## Factor1
## Religion.Important 0.837
## PubOp 0.663
## Church 0.700
## Abor 0.561
## Irreligion 0.870
## AdulteryUnac 0.378
## HomUn 0.758
## AbortionOp 0.889
## PremSexUnac 0.690
##
## Factor1
## SS loadings 4.683
## Proportion Var 0.520
##
## Test of the hypothesis that 1 factor is sufficient.
## The chi square statistic is 86.75 on 27 degrees of freedom.
## The p-value is 3.45e-08
Faith$scores
## Factor1
## Albania 0.27460352
## Argentina -0.14545910
## Armenia 1.39605320
## Australia -1.05719038
## Austria -0.85506924
## Belarus -0.15603821
## Belgium -1.40247968
## Belize 0.61591818
## Bolivia 1.13317521
## Bosnia and Herzegovina 0.69938286
## Brazil 0.74817131
## Bulgaria -0.07140885
## Canada -1.03093144
## Chile -0.04320130
## Colombia 0.87728276
## Costa Rica 0.76598928
## Croatia 0.42238665
## Cyprus 0.36849602
## Czech Republic -1.39661687
## Denmark -1.53893430
## Ecuador 1.12393797
## El Salvador 1.52338039
## Estonia -1.29150759
## Finland -1.12811617
## France -1.47303927
## Georgia 1.05818497
## Germany -1.28874922
## Greece 0.17678240
## Guatemala 1.81152952
## Honduras 0.04048214
## Hungary -0.93199542
## Iceland -1.09271411
## Ireland -1.07336030
## Italy -0.13686091
## Latvia -0.60043094
## Lithuania -0.08941334
## Mexico 0.33703252
## Moldova 0.77741428
## Montenegro 0.52593919
## Netherlands -1.20332816
## New Zealand -0.89076086
## Nicaragua 1.28617523
## North Macedonia 1.13095654
## Norway -1.42938737
## Panama 1.19877090
## Paraguay 1.44781723
## Peru 0.88068809
## Poland 0.51477570
## Portugal -0.10746666
## Romania 0.65611123
## Russia -0.13570222
## Serbia 0.20964138
## Slovakia 0.35566716
## Slovenia -0.80563264
## Spain -1.12101382
## Suriname 1.60654850
## Sweden -1.72328246
## Switzerland -0.75727603
## Ukraine 0.41079574
## United Kingdom -1.27618532
## United States -0.03895117
## Uruguay 0.25395288
## Venezuela 0.78406461
## Guyana 0.88039578
Having created two scales with sufficient reliability, let’s go back to our questions:
Below, I generate two plots. The first one pools together all countries for a scatterplot with an imposed regression line constituting the best linear fit for faith against development. We see a negative and discernible slope–in other words, that slope surely can’t be argued to be 0. I additionally provide the location of the coordinates of the United States, and a -45 degree line associated with it, from which we can identify which countries are at a comparable level of associated modernization and secularism. The answer to this question is…basically none.
The second plot creates separate regression lines for the different continents (ignoring that I’ve fancifully created a “continent” of U.S./Canada/Australia/New Zealand, all settler societies). We see that all the slopes are negative, though Latin America’s looks a bit flatter than the other ones. In other words, wealthier countries are basically more secular in general, but in Latin America, the relationship seems less strong than in other continents.
A few caveats are in order, of course:
scores <- data.frame(
Country=df$Country %>% as.vector,
Continent=df$Continent %>% as.vector,
Development=Development$scores[,1] %>% as.vector,
Faith=Faith$scores %>% as.vector
)
xintercept=scores %>% filter(Country=="United States") %>% pull(Development)
yintercept=scores %>% filter(Country=="United States") %>% pull(Faith)
p <- ggplot(scores,
aes(x=Development, y=Faith,
label=Country)) +
geom_point(aes(colour=Continent)) +
geom_text_repel() +
theme_dark() +
geom_smooth(method='lm') +
geom_vline(xintercept=scores %>% filter(Country=="United States") %>% pull(Development)) +
geom_hline(yintercept=scores %>% filter(Country=="United States") %>% pull(Faith)) +
geom_abline(intercept=xintercept+yintercept,
slope=-1)
ggplotly(p)
## `geom_smooth()` using formula 'y ~ x'
## Warning in geom2trace.default(dots[[1L]][[1L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomTextRepel() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
p <- ggplot(scores,
aes(x=Development, y=Faith,
label=Country, colour=Continent,
group=Continent)) +
geom_point() +
geom_text_repel() +
theme_dark() +
geom_smooth(method='lm') +
geom_vline(xintercept=scores %>% filter(Country=="United States") %>% pull(Development)) +
geom_hline(yintercept=scores %>% filter(Country=="United States") %>% pull(Faith)) +
geom_abline(intercept=xintercept+yintercept,
slope=-1)
ggplotly(p)
## `geom_smooth()` using formula 'y ~ x'
## Warning in geom2trace.default(dots[[1L]][[4L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomTextRepel() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
## Warning in geom2trace.default(dots[[1L]][[4L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomTextRepel() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
## Warning in geom2trace.default(dots[[1L]][[4L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomTextRepel() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
## Warning in geom2trace.default(dots[[1L]][[4L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomTextRepel() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
But wait, there’s more…
What would such a report be without a statistical model? Well, shorter, for one thing.
We’ll start with a simple model–literally, simple linear regression. To be a bit more proper, we’ll assess what effect a one “unit” change in development has on Faith. Of course, these units are constructs, not empirical quantities, so we might want to change our parlance a bit–does modernization decrease faith by a little, a lot, or not at all? Or might we see something completely unexpected? (Spoilers: No)
mod <- lm(Faith~Development, data=scores)
summary(mod)
##
## Call:
## lm(formula = Faith ~ Development, data = scores)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.17245 -0.35131 -0.00739 0.43739 0.83926
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.912e-17 6.491e-02 0.00 1
## Development -8.203e-01 6.570e-02 -12.49 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5193 on 62 degrees of freedom
## Multiple R-squared: 0.7155, Adjusted R-squared: 0.7109
## F-statistic: 155.9 on 1 and 62 DF, p-value: < 2.2e-16
plot_model(mod, type='eff')
## $Development
Survey says….surprise surprise, modernization has a negative effect on faith. A shocking, breathtaking conclusion says…nobody. The slope is certainly something to type home about. To put it in perspective, let’s examine our dependent variable:
summary(scores$Faith)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.7233 -0.9567 0.1086 0.0000 0.7688 1.8115
The minimum value is -1.7232825, and the maximum value is 1.8115295. Some fancy math tells us the range of the variable is -3.534812. The coefficient for the model tells us that going one unit up in Modernization decreases that Faith variable by -0.8203049 (thank you, R weirdness with single vs. double brackets). So, think about how much relatively you go down in faith with a single unit change in Modernization.
Ok, I’ve got a fever, and the only cure is Tylenol. No, actually, it’s…more statistical models. We saw a difference in regions in the plots above so let’s add that to the model: maybe each region has its own baseline.
mod2 <- lm(Faith~Development + Continent, data=scores)
summary(mod2)
##
## Call:
## lm(formula = Faith ~ Development + Continent, data = scores)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.18519 -0.29383 0.03051 0.45954 0.97720
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.01171 0.10863 -0.108 0.9145
## Development -0.63432 0.12590 -5.038 4.74e-06 ***
## ContinentLatin America 0.29945 0.17455 1.716 0.0915 .
## ContinentNorth America/Australasia 0.01883 0.33021 0.057 0.9547
## ContinentWestern Europe -0.33219 0.24488 -1.357 0.1801
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5097 on 59 degrees of freedom
## Multiple R-squared: 0.7391, Adjusted R-squared: 0.7214
## F-statistic: 41.79 on 4 and 59 DF, p-value: < 2.2e-16
plot_model(mod2, type='eff')
## $Development
##
## $Continent
As our plots above had suggested, Latin America seems to be a bit higher in faith than the other three regions, which are otherwise not able to be distinguished one from another. The slope of our Development variable does change a bit–some of its effect in the simple model likely owed to Latin American countries scoring lower on Development while also having lower levels of Secularization.
Are we sure e like the more complicated model? Well, beyond the R squared difference, ANOVA can do a bit for us:
anova(mod2, mod)
## Analysis of Variance Table
##
## Model 1: Faith ~ Development + Continent
## Model 2: Faith ~ Development
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 59 15.326
## 2 62 16.717 -3 -1.3911 1.7851 0.1598
One last little quibble–we may have heteroscedasticity, one of the deadly sins of ordinary least squares estimation. Just in case, since our Breusch-Pagan test is basically borderline, we’ll use some robust standard errors to make ourselves a little happier. I’m happier already, aren’t you?
bptest(mod2)
##
## studentized Breusch-Pagan test
##
## data: mod2
## BP = 6.3268, df = 4, p-value = 0.176
coeftest(mod2, vcov = sandwich::vcovHC, type = "HC0")
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.011707 0.124850 -0.0938 0.92561
## Development -0.634321 0.141371 -4.4869 3.402e-05 ***
## ContinentLatin America 0.299452 0.161438 1.8549 0.06861 .
## ContinentNorth America/Australasia 0.018826 0.305046 0.0617 0.95100
## ContinentWestern Europe -0.332195 0.268071 -1.2392 0.22018
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The results are in and we are keeping the more complicated model.
But before we go let’s try one more thing–what we did above is called a random intercept model, if we want to get a bit particular. What if we want to test random slope and random intercept? Well, we can do so fairly easily–we’ll add an interaction term for it.
mod3 <- lm(Faith~Development*Continent, data=scores)
summary(mod3)
##
## Call:
## lm(formula = Faith ~ Development * Continent, data = scores)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.03934 -0.34304 0.03864 0.29434 0.90003
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.056211 0.113302 -0.496
## Development -0.813588 0.188429 -4.318
## ContinentLatin America 0.568494 0.232673 2.443
## ContinentNorth America/Australasia 2.054041 1.977623 1.039
## ContinentWestern Europe -0.064857 0.469826 -0.138
## Development:ContinentLatin America 0.429762 0.266461 1.613
## Development:ContinentNorth America/Australasia -1.478803 1.641780 -0.901
## Development:ContinentWestern Europe -0.007565 0.412825 -0.018
## Pr(>|t|)
## (Intercept) 0.6218
## Development 6.49e-05 ***
## ContinentLatin America 0.0177 *
## ContinentNorth America/Australasia 0.3034
## ContinentWestern Europe 0.8907
## Development:ContinentLatin America 0.1124
## Development:ContinentNorth America/Australasia 0.3716
## Development:ContinentWestern Europe 0.9854
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5055 on 56 degrees of freedom
## Multiple R-squared: 0.7564, Adjusted R-squared: 0.7259
## F-statistic: 24.84 on 7 and 56 DF, p-value: 4.998e-15
plot_model(mod3)
The results are in! And, I don’t think I’m too impressed. Maybe that slope for North America/Oceania is steeper than the rest. But…maybe not? Meh. I think relative simplicity rules on this one.
Let’s turn to the bullet points.